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EXTENSION  OF  PROGRAMS  FOR  CALCULATIONS  OF  GREAT 
CIRCLE  PATHS  AND  SUNRISE-SUNSET  TIMES 

J,  H.  Crary 

A  previous  note  by  Brady  and  Crombie  [1964] 
describes  methods  and  a  program  for  calculating  and 
plotting  sunrise  and  sunset  times  at  specified 
heights  and  distances  along  a  great  circle  path 
determined  by  the  locations  of  the  end  points .   The 
program  also  plots  the  distribution  of  darkness  or 
daylight  along  the  path  versus  time. 

This  program  has  been  modified  to  remove  some 
restrictions  on  the  use  of  the  earlier  program,  and 
to  provide  the  option  of  calculating  the  short  or 
long  great  circle  path.   In  addition,  transmitter, 
receiver,  and  solar  bearings  are  also  calculated  and 
plotted.   Some  interesting  characteristics  of  the 
short  and  long  paths  from  Rugby,  England,  to  Byrd 
Station,  Antarctica,  are  shown. 

1.   Introduction 
A  previous  note  by  Brady  and  Crombie  [1964]  described  the  methods 
and  gave  the  formulas  for  the  calculation  of  sunrise  and  sunset  times 
at  specified  heights  and  distances  along  a  great  circle  path.   The 
program  previously  described  had  some  limitations  which  restricted  its 
use.   Modifications  were  made  to  the  program  to  remove  the  preceding 
restrictions  and  to  add  some  additional  features. 


These  features  are  the  calculation  of  (l)  either  the  short  or  long 
path,  (2)  the  transmitter  bearing  at  each  segment  of  the  path,  (3)  the 
solar  bearing  angle  at  sunrise  and  sunset  at  each  segment,  (k)   the  dif- 
ference angle  between  the  transmitter  and  solar  bearings,  and  (5)  the  co- 
sine of  this  difference  angle.  A  routine  was  then  added  to  plot  these 
results,  together  with  the  percentage  of  darkness  along  the  path.   The 
necessity  of  giving  the  easternmost  end  of  the  path  first  was  also  re- 
moved, since  the  choice  depends  upon  whether  the  short  or  long  path  is 
desired. 

These  modifications  were  accomplished  by  adding  three  subroutines  to 
the  original  program.  These  programs  and  subroutines  are  available  in 
Fortran  language.  A  description  of  these  will  be  given  which  includes 
the  necessary  equations  and  logical  statements  to  enable  a  person  to  cal- 
culate these  results  manually  or  by  means  of  a  computer  program,  without 
first  determining  how  to  make  the  necessary  decisions;  this  is  a  time  con- 
suming process  and  need  not  be  duplicated. 

2 .   The  Computer  Programs 
2.1.   Introduction 

The  computer  programs  consist  of  two  main  Fortran  programs,  DALITP 
and  DALITB,  together  with  the  necessary  subroutines  (PATH,  SOLANG,  DATBX, 
PLOT).   The  functions  of  these  subroutines  are:   (l)  PATH--to  calculate 
the  necessary  information  about  the  great  circle  path,  such  as  the  total 
length,  the  bearing  angles  of  the  transmitter  and  receiver,  and  the  posi- 
tion of  points  at  specified  intervals  along  the  path;  (2)  SOLANG — to 


calculate  the  solar  zenith  and  bearing  angles;  (3)  DATBX — a  function 
which  furnishes  the  date  on  which  the  calculations  are  performed  to  the 
calling  program;  and  (k)   PLOT--a  general  output  plotting  routine.  The 
functions  of  these  subroutines  are  discussed  in  2.2.  -  2.5.   The  follow- 
ing notation  will  be  used  in  the  description: 

L  =  latitude  in  degrees  (North  is  +,  South  is  -) 
0  =  longitude  in  degrees  (East  is  +,  West  is  -) 
LHA.  =  solar  local  hour  angle  in  degrees 

DECL  =  solar  declination  in  degrees 

Subscript  i  or  j  indicates  ith  or  jth  point  of  the  path 

R  indicates  receiver 

T  indicates  transmitter 

2.2.  Sunrise-Sunset  Program  (DALITP) 

This  program  is  the  modified  version  of  the  original  DALIT  program 
and  controls  the  data  card  reading  and  the  output  printing.   It  also  per- 
forms the  calculations  of  the  sunrise  and  sunset  times,  using  the  same 
equations  given  by  Brady  and  Crombie  [196U]. 

In  the  case  of  | cos  LHA|  >  1,  the  following  test  is  made: 

if  (L.  •  DECL)  <  0  the  path  is  always  dark, 
if  (Li  •  DECL)  2:  0  the  path  is  always  light. 

This  takes  care  of  both  northern  and  southern  hemisphere  conditions. 
The  data  inputs  are  the  same  as  before  except  for  the  following 


changes:   (l)  a  given  number  of  pairs  of  ionosphere  and  screen  heights, 
up  to  a  maximum  of  ten,  may  be  used,  and  (2)  the  option  of  calculating 
the  long  or  short  path  is  provided. 

2.3.   Sunrise -Sunset  Program  with  Bearing  Calculations  (DALITB) 

This  program  is  the  same  as  DALITP,  with  the  addition  of  the  sub- 
routine SOLANG  to  calculate  the  solar  bearing  angles  between  the  solar 
and  transmitter  bearing,  and  the  cosine  of  this  difference.  The  PLOT 
subroutine  then  is  used  to  plot  the  solar  bearings,  the  cosine  of  the 
difference  angle,  and  the  distance  along  the  path  where  the  terminator* 
crosses,  as  a  function  of  time. 

2.1+.   Great  Circle  Path  Subroutine  (PATH) 

This  subroutine  is  called  by  the  main  program  to  do  the  great  circle 
calculations.  The  program  has  several  options  to  calculate  for  the  short 
or  long  path:  (l)  distance  between  two  specific  points;  (2)  latitude, 
longitude,  and  bearing  of  transmitter  at  given  intervals  on  a  given  path; 
(3)  coordinates  of  reflection  points  for  an  n-ref lection  ray  between  given 
points;  or  (k)  coordinates  of  the  end  points  of  a  path  when  given  one  end 
and  the  jth  point  on  an  n-ref lection  path. 

The  option  of  being  able  to  calculate  either  the  short  or  the  long 
great  circle  path  between  two  specified  points  requires  several  logical 


*  Terminator--a  useful  astronomical  term  for  the  line  separating  the  dark 
and  sunlit  portions  of  a  planet,  applied  here  at  the  ionospheric  height 
specified. 


decisions  to  be  made  at  different  points  during  the  calculation  in  order 
to  insure  that  the  proper  solution  is  being  obtained.   The  spherical 
triangles  involved  in  these  solutions  are  illustrated  in  figure  1.   The 
polar  angle  of  the  spherical  triangle  is  the  difference  in  longitude  be- 
tween the  two  points  on  the  path  being  considered.  The  algebraic  sign  of 
the  difference  is  not  used  in  the  calculations.   Four  cases  can  actually 
be  defined,  depending  upon  whether  the  short  or  long  path  is  desired,  and 
whether  the  transmitter  or  starting  end  of  the  path  is  the  eastern-most  or 
western-most  point  of  the  path.  This  decision  depends  upon  whether  the 
short  or  long  path  is  desired. 

In  addition,  decisions  must  be  made  as  to  which  of  the  possible 
solutions  to  equations  involving  cos"1  will  yield  the  desired  path.   It 
is  also  necessary  to  account  for  some  of  the  cases  of  the  degenerate 
spherical  triangle,  such  as  where  both  points  of  the  path  lie  on  a 
meridian,  etc.  For  this  reason,  the  necessary  equations  and  logic  state- 
ments are  included,  permitting  a  straightforward  calculation  of  the 
desired  path.  The  flagl  mentioned  in  the  following  description  is  a 
logical  variable  which  indicates  a  decision  made  on  the  basis  of  the 
orientation  of  transmitter  and  receiver  and  whether  the  short  or  long 
path  is  desired. 

(l)        Calculate  the  polar  angle  or  longitude  difference,  |^  -  $~\ 

(see  figure  l) . 
(2.1)       If  |fzL  -  <f>   |  <  180°,  continue;  if  >  l80°  ,  go  to  (2.2). 

(2.1.1)   If  short  path,  <f>  =   360  -  \<f>     -  </>    |  and  set  flagl  =  2,  go  to 
(3). 


(2.1.2)  If  long  path,  0  =  |0T  -  0R|  ,  and  set  flagl  =  1,  go  to  (3). 

(2.2)  If  |0T  -  0R|  >  180°. 

(2.2.1)  If  short  path,  0  =  |0T  -  0R| ,  and  set  flagl  =  1,  go  to  (3). 

(2.2.2)  If  long  path,  0  =  36O  -  |0T  -  0R|  ,  and  set  flagl  =  2. 

(3)         Calculate  D,  the  total  path  distance. 

(3.1)  If  cos  LR  ±   0,  go  to  (3.2);  if  =  0,  continue.   If  cos  %  =  0, 
then  DRAD  =  |LT  -  LJ  *  tt/180  radians. 

(3.1.1)  If  LR  =  0,  impossible  branch,  error  message,  exit. 

(3.1.2)  If  LR  <  0,  T  =  tt  rad,  where  T  is  the  angle  of  the  spherical 
path  triangle  at  the  transmitter. 

(3.1.3)  If  LR  >  0,  T  =  0. 

(3.1.4)  If  long  path  DRAD  -  DRAD  -  2tt.  D  =  DRAD  *  R  where  R  is  the 
effective  earth  radius.  R  is  taken  as  the  temperate  latitude 
average  of  the  radii  determined  from  the  polar  and  equatorial 
circumferences,  or  (3  Rea  +  R^ol)/^*  r^ae   calculation  then 
proceeds  to  (k) . 

(3.2)  If  cos  LR  ^  0,  DRAD  =  cos"1  [sin  Lm  sin  L-^  +  cos  Lrp  cos  LR 
cos  0]  from  the  cosine  law  for  the  spherical  triangle  shown 
in  figure  1. 

(3.2.1)  If  sin  (DRAD)  =  0,  then  end  points  are  antipodal,  and  there 
are  an  infinite  number  of  solutions.  This  causes  a  message 
to  be  printed  and  the  calculation  is  terminated. 

(3.2.2)  If  sin  (DRAD)  ^  0:   for  long  path  DRAD  =  2tt  -  DRAD;  for  both 
paths  D  =  DRAD  •  R.   If  only  distance  is  desired,  return  to 
main  program. 


(3.2.3)  Determine  number  of  segments.   Complete  calculation  of  angle 
T. 

(4.1)  If  cos  LT  =  0,  continue;  if  #   0,  go  to  (4.3) • 
(U.l.l)   If  L  =  0,  error  message,  exit. 

(4.1.2)  If  LT  <  0,  T  =  0,  go  to  (4.1.4). 

(4.1.3)  If  LT  >  0,  T  =  tt. 

(4.1.4)  Then  Lj  =  LT  -  (LT/|LT|)  D-/R,    B  =  90  -  (LT/|LT| )  90,  0  = 

0R;  go  "to  (4.3). 

(4.2)  If  cos  LT  f   0. 

(4.3)  T  =  cos"1  [(sin  L-n  -  cos  D  sin  L  )/(cos  L<p  sin  D)]  from  the 
law  of  cosines  again  applied  to  the  path  triangle  of  figure  1. 

(5)  Calculate  the  latitudes  at  the  jth  segment  along  the  p~i,th: 

L .  =  90   -  cos"1 (sin  Lm  cos  D^  +  cos  L  sin  D.  cos  T) .   This 
is  another  variation  of  the  law  of  cosines  applied  to  the  tri- 
angle of  the  pole,  point  j  and  the  transmitter.  The  use  of 
90  -  cos"1  instead  of  sin"1  removes  an  ambiguity  in  the  result, 

(5.1)  If  cos  L.  =0,  A0  =  0,  since  point  j  is  at  a  pole.   Here  A0 

J 

is  the  difference  in  longitude  between  the  transmitter  and 
point  j.   Go  to  (6.3) . 

(5.2)  If  cos  L.  /   0. 

J 

(6)  Calculate  the  longitudes. 

(6.l)       Cos  A0  =  (cos  D.  -  sin  L^  sin  L.)/(cos  L™  cos  L.).   This  also 
comes  from  the  law  of  cosines. 


(6.2) 


Now,  using  the  law  of  sines  for  sin  A0:   sin  A0  =  (sin  T  sin 


Dp/Ccos  Lj). 

(6.2.1)  If  sin  A0  ^  0,  A0  =  cos"1 (cos  A0),  go  to  (6.3). 

(6.2.2)  If  sin  A0  <  0,  A0  =  36O  -  cos"1 (cos  A0). 

(6.3)  If  (LT  -  L  )  £  0,  continue;  if  <  0,  go  to  (6.4). 

(6.3.1)  If  flagl  =  1  [see  (2.l)-(2  .2)] ,  0j  =  0T  -  A0,  go  to  (7). 

(6.3.2)  If  flagl  =  2,  0  =  0T  +  A0,  go  to  (7). 
(6.1+)  If  (LT  -  LR)  <  0  and 

(6.4.1)  If  flagl  =  1  use  (6.3.2),  go  to  (7). 

(6.4.2)  If  flagl  =  2  use  (6.3.1). 


(7) 
(7.1) 


Check  longitude  for  proper  hemisphere 

If  |0.|  -  180  £  0,  continue;  if  >  0,  go  to  (8). 

J 


(7.1.1)   If  0-    <    0,  0.    =    0.    +  360°,  gO  to  (8). 


(7.2) 


If  |0,j  ^  0,  0   =  360°  +  0.. 


(8). 

(8.1) 

(8.2) 

(8.3) 


Find  transmitter -bearing  B.. 


If  cos  \  =  t,   B  =  90  +  (I|R)/C1LR!0  •  90,  go  to  (9). 

Or  if  cos  L_  ^   0,  B.  =  cos-1  [(sin  Lm  -  sin  L.  cos  D.)/ 

(cos  L.  sin  D.)]  again  from  the  law  of  cosines. 
3              3 

If  (^  -  <f>  )   <  0,  continue;  if  >  0,  go  to  (8.4). 


(8.3.1)  Flagl  =  1  [see  (2.l)-(2.2)]  B.  =  36O  -  B.,  go  to  (9) 

J  J 

(8.3.2)  Or  if  flagl  =2,  B.  =  B.,  go  to  (9). 

J        J 

(8.4)      if  (^  -  ^)  *  0. 

(8.4.1)  And  flagl  =  1,  use  (8.3.2),  go  to  (9). 

(8.4.2)  Or  if  flag  1=2,  use  (8.3.1). 
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(9)  Find  receiver  "bearing  from  transmitter. 

(9.1)  If  (0T  -  0R)  <  0,  continue;  if  S  0,  go  to  (9.2). 

(9.1.1)  And  flagl  =  1,  BR  =  T,  go  to  (10). 

(9.1.2)  Or  if  flagl  =2,  BR=  36O  -  T,  go  to  (10). 

(9.2)  If  (0T  -  0R)  £  0. 

(9.2.1)  And  flagl  =  1  use  (9.1.2),  go  to  (10). 

(9.2.2)  Or  if  flagl  =  2  use  (9.I.1). 

(10)  Since  the  transmitter  "bearing  is  undefined  at  the  transmitter 
we  use: 

(10.1)  Bj.  =  Bg  +  180. 

(10.2)  If  (B1  -  360)  £  0,  B-l  =  B1  -  360,  go  to  (ll), 

(10.3)  Or  if  (Bx  -  360)  <  0,  B1  =  Bj_. 

(11)  The  subroutine  returns  to  the  main  program  here. 

2.5.  Solar  Angle  Subroutine  (SOLANG) 


This  subroutine  finds  the  solar  bearing  angles  S.  and  the  zenith 
angles  \±. 

The  inputs  are  an  array  of  N  sets  of  path  latitude  (XLAT)  and  longi- 
tude (XLON) ,  solar  declination  (DECL) ,  and  Greenwich  hour  angle  (GHA), 
where  N  ^  1000.   The  variable  NTIMES  may  be  used  to  specify  whether  there 
is  only  one  pair  of  values  (NTIMES  =  l)  for  GHA  and  DECL,  or  whether  there 
is  a  pair  for  each  point  (NTIMES  >  l). 

The  logic  statements  given  here  are  used  to  find  the  proper  solution 
to  equations  involving  the  arc  cosine  function.   It  is  necessary  to  pick 


the  solution  corresponding  to  the  desired  triangle.   It  is  also  necessary 
to  make  decisions  so  that  the  solar  bearing  angle  will  follow  the  conven- 
tion of  being  measured  clockwise  from  North  and  the  solar  zenith  angle 
will  lie  between  -l80°  and  +l80° . 

The  calculation  then  proceeds  as  follows,  using  the  solar  triangle 
of  figure  2 . 

(1)  Find  the  solar  local  hour  angle  SLHA  =  GHA  +  $.. 

(1.1)  If  SEHA  <  0,  then  SLHA  =  SLHA  +  36O,  go  to  (l.3). 

(1.2)  If  SLHA  >  0  then 

(1.3)  If   (SLHA  -   360)   >  0,    continue;    if  <  0,   go  to   (1.4). 
(1.3.1)  SLHA  =  SLHA  -   36O,   go  to    (2). 

(1.4)  If    (SLHA  -   360)   <  0,    SLHA  =  SLHA. 

(2)  Calculate  v.    =  cos-1 (sin  L.    sin  DECL  +  cos  L.    cos  DECL 

M.  1  1 

cos  SLHA)  from  the  law  of  cosines.   Next  find  the  bearing 
angle. 

(2.1)  If  I L.|  -  90  =  0,  continue;  if  ^   0,  gc  to  (2.2) . 
(2.1.1)   S.  =  SLHA,  go  to  (2.5). 

(2.2)  If  |L.|  -  90  ^  0. 

(2.3)  If  (xlxl  "  l8°)  =  °>  continue;  if  ^  0,  go  to  (2.4). 

(2.3.1)  And  if  L.  -  DECL  <  0,  S.  =  l80° ,  go  to  (2.5). 

(2.3.2)  Or  if  L.  -  DECL  >  0,  S.  =  0°,  go  to  (2.5). 

(2.4)  If  x(|x|  -  180)  +   0. 

(2.4.1)   S.  =  cos-1 [(sin  DECL  -  sin  L.  cos  x)/(cos  L.  sin  y)]   from  the 
law  of  cosines. 

(2.5)  If  (SLHA  -  180)  <  0,  S.  =  360  -  S.,  go  to  (3). 
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(2.6)       Or  if  (SLHA  -  180)  ^  0,  S^^  =  S±. 
(3)        Return  to  the  main  program. 

3.  Some  Illustrations  of  the  Results  for  Long  and  Short  Paths 

Recordings  of  VLF  signals  propagating  over  long  VLF  paths  (>  15,000 
km  or  so)  often  show  diurnal  variations  in  phase  and  amplitude  which  are 
indicative  of  at  least  some  energy  coming  over  the  long  great  circle  path 
[Crombie,  1958] ,  particularly  when  the  short  path  is  completely  or  almost 
completely  sunlit.  This  makes  the  variation  of  daylight  or  darkness  over 
both  the  long  and  short  paths  especially  interesting. 

Figures  3  and  k   show  typical  sets  of  the  computer  output  for  sunrise 
and  sunset  calculations  for  the  short  and  long  path  from  GBR,  Rugby,  Eng- 
land, to  Byrd  Station,  Antarctica,  for  April  15,  1962 .  The  upper  set  in 
each  figure  is  taken  for  an  ionospheric  height  of  90  km,  and  a  screen 
height  of  27  km,  or  x  ~  98^ •   The  lower  set  in  each  figure  is  for  ground 
level  or  x  =  90  .   These  figures  illustrate  several  interesting  points. 
First  is  the  essentially  complementary  nature  of  the  two  paths,  i.e.,  one 
is  becoming  light  when  the  other  is  becoming  dark,  and  vice  versa.  An<r 
other  is  the  extremely  rapid  transition, especially  at  ground  level, around 
1900  UT.   This  could  be  expected  to  cause  a  very  rapid  transition  of  the 
predominant  component  of  the  received  signal  from  the  long  path  to  the 
short  path.  This  effect  occurs  when  the  great  circle  joining  the  trans- 
mitter and  receiver  nearly  coincides  with  the  terminator.   It  is  indicated 
by  a  solar  and  transmitter  bearing  difference  close  to  90°  or  270°  at  sun- 
rise or  sunset. 
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These  effects  are  also  illustrated  in  figures  5  and  6.   These  show 
the  output  plot  of  percentage  of  path  in  darkness,  distance  along  path  of 
terminator  location,  the  solar  bearing  angle,  and  the  cosine  of  the  solar- 
transmitter  bearing  angle  difference  at  sunrise  and  sunset.   The  rapid 
transition  occurs  when  the  bearing  angle  at  which  the  path  crosses  the 
equator  (9.63° )  is  nearly  equal  in  magnitude  to  the  solar  declination 
(9.6^°).   The  bearing  at  which  the  path  crosses  the  equator  is  also  equal 
to  the  complement  of  the  maximum  latitude  (80.37°)  on  the  great  circle 
path.   The  bearing  difference  angle  is  interesting,  not  only  because  of 
its  influence  on  the  time  interval  during  which  the  night-day  transition 
is  made,  but  also  because  experimental  data  indicate  that  the  angle  at 
which  the  terminator  crosses  the  path  appears  to  influence  the  amount  of 
mode  conversion  at  the  discontinuity  in  the  waveguide  height  caused  by 
the  terminator  [Crombie,  I96U] . 

Figure  6  illustrates  an  interesting  characteristic  of  some  paths  at 
certain  times .  Note  that  the  path  at  ground  level  remains  about  80$  dark 
until  I85O  UT.  The  percent  of  darkness  then  drops  to  18$  in  about  fifteen 
minutes.   The  percentage  in  darkness  at  the  ionospheric  height  changes 
gradually  from  about  1^30  UT.   It  begins  to  drop  rapidly  after  1700  and 
reaches  zero  at  1825  UT.   It  then  rises  from  zero  to  10$  between  1950  and 
2120  UT. 

A  corresponding  phase  retardation  at  about  the  proper  time  (1958  UT 
in  this  case)  on  a  path  showing  such  behavior  may  enable  one  to  determine 
the  critical  value  of  x,  at  which  the  ionizing  radiation  first  affects 
the  D  layer. 
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Figure 


Geometry  of  a  Great  Circle  Path 
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The  relationship  of  the  solar  local  hour  angle,  zenith  angle,  declination  and  bearing  angle, 
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Figure  3.  Light  and  Darkness  Plot  vs.  UT  from  the  Short  Path  from  GBR,  Rugby,  England,  to 
Byrd  Station,  Antarctica  for  15  April  (1962),  for  (a)  ground  or  y  =  90  ,  ang 
(b)  ionospheric  reflection  height .ss  90  km,  screen  height  =  27  km,  or  \  =  9&   . 
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Figure  4.   Light  and  Darkness  Plots  vs.  UT  for  Long  Path  from  GBR,  Rugby,  Engl 
Byrd  Station,  Antarctica  for  15  April  (1962),  for  (a)  ground  or  x 
(b)  ionospheric  reflection  height  -  90  km,  screen  height  =  27  km,  o 
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